Picture1.png

Picture1.png

Exploration¶

Load Images¶

In [20]:
#import imageio.v2 as imageio
import imageio
import pydicom
from pydicom.data import get_testdata_file
import numpy as np

# load a single DICOM image from the scan volume

ct_file = get_testdata_file("CT_small.dcm")
mr_file = get_testdata_file("MR_small.dcm")
im_ct = imageio.imread(ct_file, progress=True) 
im_mr = imageio.imread(mr_file, progress=True)
C:\Users\nasta\AppData\Local\Temp/ipykernel_25900/3685182357.py:11: DeprecationWarning: Starting with ImageIO v3 the behavior of this function will switch to that of iio.v3.imread. To keep the current behavior (and make this warning disappear) use `import imageio.v2 as imageio` or call `imageio.v2.imread` directly.
  im_ct = imageio.imread(ct_file, progress=True)
C:\Users\nasta\AppData\Local\Temp/ipykernel_25900/3685182357.py:12: DeprecationWarning: Starting with ImageIO v3 the behavior of this function will switch to that of iio.v3.imread. To keep the current behavior (and make this warning disappear) use `import imageio.v2 as imageio` or call `imageio.v2.imread` directly.
  im_mr = imageio.imread(mr_file, progress=True)

Medical Image Metadata¶

Comprising : Patient Demographics and Acquisition Information

In [6]:
ds_ct = pydicom.dcmread(ct_file)
ds_ct
Out[6]:
Dataset.file_meta -------------------------------
(0002, 0000) File Meta Information Group Length  UL: 192
(0002, 0001) File Meta Information Version       OB: b'\x00\x01'
(0002, 0002) Media Storage SOP Class UID         UI: CT Image Storage
(0002, 0003) Media Storage SOP Instance UID      UI: 1.3.6.1.4.1.5962.1.1.1.1.1.20040119072730.12322
(0002, 0010) Transfer Syntax UID                 UI: Explicit VR Little Endian
(0002, 0012) Implementation Class UID            UI: 1.3.6.1.4.1.5962.2
(0002, 0013) Implementation Version Name         SH: 'DCTOOL100'
(0002, 0016) Source Application Entity Title     AE: 'CLUNIE1'
-------------------------------------------------
(0008, 0005) Specific Character Set              CS: 'ISO_IR 100'
(0008, 0008) Image Type                          CS: ['ORIGINAL', 'PRIMARY', 'AXIAL']
(0008, 0012) Instance Creation Date              DA: '20040119'
(0008, 0013) Instance Creation Time              TM: '072731'
(0008, 0014) Instance Creator UID                UI: 1.3.6.1.4.1.5962.3
(0008, 0016) SOP Class UID                       UI: CT Image Storage
(0008, 0018) SOP Instance UID                    UI: 1.3.6.1.4.1.5962.1.1.1.1.1.20040119072730.12322
(0008, 0020) Study Date                          DA: '20040119'
(0008, 0021) Series Date                         DA: '19970430'
(0008, 0022) Acquisition Date                    DA: '19970430'
(0008, 0023) Content Date                        DA: '19970430'
(0008, 0030) Study Time                          TM: '072730'
(0008, 0031) Series Time                         TM: '112749'
(0008, 0032) Acquisition Time                    TM: '112936'
(0008, 0033) Content Time                        TM: '113008'
(0008, 0050) Accession Number                    SH: ''
(0008, 0060) Modality                            CS: 'CT'
(0008, 0070) Manufacturer                        LO: 'GE MEDICAL SYSTEMS'
(0008, 0080) Institution Name                    LO: 'JFK IMAGING CENTER'
(0008, 0090) Referring Physician's Name          PN: ''
(0008, 0201) Timezone Offset From UTC            SH: '-0500'
(0008, 1010) Station Name                        SH: 'CT01_OC0'
(0008, 1030) Study Description                   LO: 'e+1'
(0008, 1090) Manufacturer's Model Name           LO: 'RHAPSODE'
(0009, 0010) Private Creator                     LO: 'GEMS_IDEN_01'
(0009, 1001) [Full fidelity]                     LO: 'GE_GENESIS_FF'
(0009, 1002) [Suite id]                          SH: 'CT01'
(0009, 1004) [Product id]                        SH: 'HiSpeed CT/i'
(0009, 1027) [Image actual date]                 SL: 862399669
(0009, 1030) [Service id]                        SH: ''
(0009, 1031) [Mobile location number]            SH: ''
(0009, 10e6) [Genesis Version - now]             SH: '05'
(0009, 10e7) [Exam Record checksum]              UL: 973283917
(0009, 10e9) [Actual series data time stamp]     SL: 862399669
(0010, 0010) Patient's Name                      PN: 'CompressedSamples^CT1'
(0010, 0020) Patient ID                          LO: '1CT1'
(0010, 0030) Patient's Birth Date                DA: ''
(0010, 0040) Patient's Sex                       CS: 'O'
(0010, 1002)  Other Patient IDs Sequence  2 item(s) ---- 
   (0010, 0020) Patient ID                          LO: 'ABCD1234'
   (0010, 0022) Type of Patient ID                  CS: 'TEXT'
   ---------
   (0010, 0020) Patient ID                          LO: '1234ABCD'
   (0010, 0022) Type of Patient ID                  CS: 'TEXT'
   ---------
(0010, 1010) Patient's Age                       AS: '000Y'
(0010, 1030) Patient's Weight                    DS: '0.0'
(0010, 21b0) Additional Patient History          LT: ''
(0011, 0010) Private Creator                     LO: 'GEMS_PATI_01'
(0011, 1010) [Patient Status]                    SS: 0
(0018, 0010) Contrast/Bolus Agent                LO: 'ISOVUE300/100'
(0018, 0022) Scan Options                        CS: 'HELICAL MODE'
(0018, 0050) Slice Thickness                     DS: '5.0'
(0018, 0060) KVP                                 DS: '120.0'
(0018, 0088) Spacing Between Slices              DS: '5.0'
(0018, 0090) Data Collection Diameter            DS: '480.0'
(0018, 1020) Software Versions                   LO: '05'
(0018, 1040) Contrast/Bolus Route                LO: 'IV'
(0018, 1100) Reconstruction Diameter             DS: '338.6716'
(0018, 1110) Distance Source to Detector         DS: '1099.3100585938'
(0018, 1111) Distance Source to Patient          DS: '630.0'
(0018, 1120) Gantry/Detector Tilt                DS: '0.0'
(0018, 1130) Table Height                        DS: '133.699997'
(0018, 1150) Exposure Time                       IS: '1601'
(0018, 1151) X-Ray Tube Current                  IS: '170'
(0018, 1152) Exposure                            IS: '170'
(0018, 1160) Filter Type                         SH: 'LARGE BOWTIE FIL'
(0018, 1190) Focal Spot(s)                       DS: '0.7'
(0018, 1210) Convolution Kernel                  SH: 'STANDARD'
(0018, 5100) Patient Position                    CS: 'FFS'
(0019, 0010) Private Creator                     LO: 'GEMS_ACQU_01'
(0019, 1002) [Detector Channel]                  SL: 912
(0019, 1003) [Cell number at Theta]              DS: '373.75'
(0019, 1004) [Cell spacing]                      DS: '1.0166'
(0019, 100f) [Horiz. Frame of ref.]              DS: '955.799988'
(0019, 1011) [Series contrast]                   SS: 2
(0019, 1013) [Start number for baseline]         SS: 0
(0019, 1014) [End number for baseline]           SS: 0
(0019, 1015) [Start number for enhanced scans]   SS: 0
(0019, 1016) [End number for enhanced scans]     SS: 0
(0019, 1017) [Series plane]                      SS: 2
(0019, 1018) [First scan ras]                    LO: 'S'
(0019, 1019) [First scan location]               DS: '7.79187'
(0019, 101a) [Last scan ras]                     LO: 'I'
(0019, 101b) [Last scan loc]                     DS: '-320.197968'
(0019, 101e) [Display field of view]             DS: '0.0'
(0019, 1023) [Table Speed [mm/rotation]]         DS: '5.0'
(0019, 1024) [Mid Scan Time [sec]]               DS: '17.784578'
(0019, 1025) [Mid scan flag]                     SS: 1
(0019, 1026) [Tube Azimuth [degree]]             SL: 0
(0019, 1027) [Rotation Speed [msec]]             DS: '1.0'
(0019, 102a) [x-ray On position]                 DS: '178.079926'
(0019, 102b) [x-ray Off position]                DS: '3994.299316'
(0019, 102c) [Number of triggers]                SL: 10431
(0019, 102e) [Angle of first view]               DS: '-718.079956'
(0019, 102f) [Trigger frequency]                 DS: '984.0'
(0019, 1039) [SFOV Type]                         SS: 16
(0019, 1040) [Stat recon flag]                   SS: 0
(0019, 1041) [Compute type]                      SS: 1
(0019, 1042) [Segment Number]                    SS: 0
(0019, 1043) [Total Segments Required]           SS: 0
(0019, 1044) [Interscan delay]                   DS: '1.0'
(0019, 1047) [View compression factor]           SS: 1
(0019, 104a) [Total no. of ref channels]         SS: 6
(0019, 104b) [Data size for scan data]           SL: 10431
(0019, 1052) [Recon post proc. Flag]             SS: 1
(0019, 1057) [CT water number]                   SS: -95
(0019, 1058) [CT bone number]                    SS: 0
(0019, 105e) [Number of channels (1...512)]      SL: 763
(0019, 105f) [Increment between channels]        SL: 1
(0019, 1060) [Starting view]                     SL: 1969
(0019, 1061) [Number of views]                   SL: 1576
(0019, 1062) [Increment between views]           SL: 1
(0019, 106a) [Dependent on #views processed]     SS: 4
(0019, 106b) [Field of view in detector cells]   SS: 852
(0019, 1070) [Value of back projection button]   SS: 1
(0019, 1071) [Set if fatq estimates were used]   SS: 0
(0019, 1072) [Z chan avg over views]             DS: '0.0'
(0019, 1073) [Avg of left ref chans over views]  DS: '0.0'
(0019, 1074) [Max left chan over views]          DS: '0.0'
(0019, 1075) [Avg of right ref chans over views] DS: '0.0'
(0019, 1076) [Max right chan over views]         DS: '0.0'
(0019, 10da) [Reference channel used]            SS: 0
(0019, 10db) [Back projector coefficient]        DS: '0.0'
(0019, 10dc) [Primary speed correction used]     SS: 1
(0019, 10dd) [Overrange correction used]         SS: 1
(0019, 10de) [Dynamic Z alpha value]             DS: '0.0'
(0020, 000d) Study Instance UID                  UI: 1.3.6.1.4.1.5962.1.2.1.20040119072730.12322
(0020, 000e) Series Instance UID                 UI: 1.3.6.1.4.1.5962.1.3.1.1.20040119072730.12322
(0020, 0010) Study ID                            SH: '1CT1'
(0020, 0011) Series Number                       IS: '1'
(0020, 0012) Acquisition Number                  IS: '2'
(0020, 0013) Instance Number                     IS: '1'
(0020, 0032) Image Position (Patient)            DS: [-158.135803, -179.035797, -75.699997]
(0020, 0037) Image Orientation (Patient)         DS: [1.000000, 0.000000, 0.000000, 0.000000, 1.000000, 0.000000]
(0020, 0052) Frame of Reference UID              UI: 1.3.6.1.4.1.5962.1.4.1.1.20040119072730.12322
(0020, 0060) Laterality                          CS: ''
(0020, 1040) Position Reference Indicator        LO: 'SN'
(0020, 1041) Slice Location                      DS: '-77.2040634155'
(0020, 4000) Image Comments                      LT: 'Uncompressed'
(0021, 0010) Private Creator                     LO: 'GEMS_RELA_01'
(0021, 1003) [Series from which Prescribed]      SS: 2
(0021, 1005) [Genesis Version - now]             SH: '05'
(0021, 1007) [Series Record checksum]            UL: 1605775145
(0021, 1015) [Unknown]                           US: 24078
(0021, 1016) [Unknown]                           SS: 2
(0021, 1018) [Genesis version - Now]             SH: '05'
(0021, 1019) [Acq recon record checksum]         UL: 750675506
(0021, 1037) [Screen Format]                     SS: 16
(0021, 104a) [Anatomical reference for scout]    LO: ''
(0021, 1090) [Tube focal spot position]          SS: 7400
(0021, 1091) [Biopsy position]                   SS: 0
(0021, 1092) [Biopsy T location]                 FL: 0.0
(0021, 1093) [Biopsy ref location]               FL: 0.0
(0023, 0010) Private Creator                     LO: 'GEMS_STDY_01'
(0023, 1070) [Start time(secs) in first axial]   FD: 862399761.111079
(0023, 1074) [No. of updates to header]          SL: 1
(0023, 107d) [Indicates study has complete info  SS: 0
(0025, 0010) Private Creator                     LO: 'GEMS_SERS_01'
(0025, 1006) [Last pulse sequence used]          SS: 0
(0025, 1007) [Images in Series]                  SL: 44
(0025, 1010) [Landmark Counter]                  SL: 0
(0025, 1011) [Number of Acquisitions]            SS: 0
(0025, 1017) [Series Complete Flag]              SL: 0
(0025, 1018) [Number of images archived]         SL: 0
(0025, 1019) [Last image number used]            SL: 4
(0025, 101a) [Primary Receiver Suite and Host]   SH: ''
(0027, 0010) Private Creator                     LO: 'GEMS_IMAG_01'
(0027, 1006) [Image archive flag]                SL: 1
(0027, 1010) [Scout Type]                        SS: 0
(0027, 101c) [Vma mamp]                          SL: 150
(0027, 101d) [Vma phase]                         SS: 1
(0027, 101e) [Vma mod]                           SL: 24
(0027, 101f) [Vma clip]                          SL: 129
(0027, 1020) [Smart scan ON/OFF flag]            SS: 1
(0027, 1030) [Foreign Image Revision]            SH: ''
(0027, 1035) [Plane Type]                        SS: 2
(0027, 1040) [RAS letter of image location]      SH: 'I'
(0027, 1041) [Image location]                    FL: -77.20406341552734
(0027, 1042) [Center R coord of plane image]     FL: -11.199999809265137
(0027, 1043) [Center A coord of plane image]     FL: 9.699999809265137
(0027, 1044) [Center S coord of plane image]     FL: -75.69999694824219
(0027, 1045) [Normal R coord]                    FL: 0.0
(0027, 1046) [Normal A coord]                    FL: 0.0
(0027, 1047) [Normal S coord]                    FL: -1.0
(0027, 1048) [R Coord of Top Right Corner]       FL: -180.53579711914062
(0027, 1049) [A Coord of Top Right Corner]       FL: 179.03579711914062
(0027, 104a) [S Coord of Top Right Corner]       FL: -75.69999694824219
(0027, 104b) [R Coord of Bottom Right Corner]    FL: -180.53579711914062
(0027, 104c) [A Coord of Bottom Right Corner]    FL: -159.63580322265625
(0027, 104d) [S Coord of Bottom Right Corner]    FL: -75.69999694824219
(0027, 1050) [Scan Start Location]               FL: -63.19999694824219
(0027, 1051) [Scan End Location]                 FL: -116.20304870605469
(0027, 1052) [RAS letter for side of image]      SH: 'L'
(0027, 1053) [RAS letter for anterior/posterior] SH: 'A'
(0027, 1054) [RAS letter for scout start loc]    SH: 'I'
(0027, 1055) [RAS letter for scout end loc]      SH: 'I'
(0028, 0002) Samples per Pixel                   US: 1
(0028, 0004) Photometric Interpretation          CS: 'MONOCHROME2'
(0028, 0010) Rows                                US: 128
(0028, 0011) Columns                             US: 128
(0028, 0030) Pixel Spacing                       DS: [0.661468, 0.661468]
(0028, 0100) Bits Allocated                      US: 16
(0028, 0101) Bits Stored                         US: 16
(0028, 0102) High Bit                            US: 15
(0028, 0103) Pixel Representation                US: 1
(0028, 0120) Pixel Padding Value                 SS: -2000
(0028, 1052) Rescale Intercept                   DS: '-1024.0'
(0028, 1053) Rescale Slope                       DS: '1.0'
(0029, 0010) Private Creator                     LO: 'GEMS_IMPS_01'
(0029, 1004) [Lower range of Pixels1]            SL: 0
(0029, 1005) [Lower range of Pixels1]            DS: '0.0'
(0029, 1006) [Lower range of Pixels1]            DS: '0.0'
(0029, 1007) [Lower range of Pixels1]            SL: 87
(0029, 1008) [Lower range of Pixels1]            SH: ''
(0029, 1009) [Lower range of Pixels1]            SH: ''
(0029, 100a) [Lower range of Pixels1]            SS: 764
(0029, 1026) [Version of the hdr struct]         SS: 2
(0029, 1034) [Advantage comp. Overflow]          SL: 0
(0029, 1035) [Advantage comp. Underflow]         SL: 0
(0043, 0010) Private Creator                     LO: 'GEMS_PARM_01'
(0043, 1010) [Window value]                      US: 400
(0043, 1011) [Total input views]                 US: 10431
(0043, 1012) [X-ray chain]                       SS: [14, 2, 3]
(0043, 1013) [Decon kernel parameters]           SS: [107, 21, 4, 2, 20]
(0043, 1014) [Calibration parameters]            SS: [4, 4, 5]
(0043, 1015) [Total output views]                SS: 10431
(0043, 1016) [Number of overranges]              SS: 0
(0043, 1017) [IBH image scale factors]           DS: '0.095'
(0043, 1018) [BBH coefficients]                  DS: [0.085000, 1.102000, 0.095000]
(0043, 1019) [Number of BBH chains to blend]     SS: 350
(0043, 101a) [Starting channel number]           SL: 7
(0043, 101b) [Ppscan parameters]                 SS: [0, 0, 0, 0, 0]
(0043, 101c) [GE image integrity]                SS: 0
(0043, 101d) [Level value]                       SS: 40
(0043, 101e) [Delta Start Time [msec]]           DS: '2.0'
(0043, 101f) [Max overranges in a view]          SL: 0
(0043, 1020) [Avg overranges all views]          DS: '0.0'
(0043, 1021) [Corrected after glow terms]        SS: 0
(0043, 1025) [Reference channels]                SS: [1, 2, 3, 748, 749, 750]
(0043, 1026) [No views ref chans blocked]        US: [0, 1, 1, 0, 0, 0]
(0043, 1027) [Scan Pitch Ratio]                  SH: '/1.0:1'
(0043, 1028) [Unique image iden]                 OB: Array of 80 elements
(0043, 1029) [Histogram tables]                  OB: Array of 2068 elements
(0043, 102a) [User defined data]                 OB: Array of 40 elements
(0043, 102b) [Private Scan Options]              SS: [4, 4, 0, 0]
(0043, 1031) [Recon Center Coordinates]          DS: [-11.200000, 9.700000]
(0043, 1040) [Trigger on position]               FL: 178.07992553710938
(0043, 1041) [Degree of rotation]                FL: 3816.219482421875
(0043, 1042) [DAS trigger source]                SL: 0
(0043, 1043) [DAS fpa gain]                      SL: 0
(0043, 1044) [DAS output source]                 SL: 1
(0043, 1045) [DAS ad input]                      SL: 1
(0043, 1046) [DAS cal mode]                      SL: 3
(0043, 1047) [DAS cal frequency]                 SL: -1
(0043, 1048) [DAS reg xm]                        SL: 1
(0043, 1049) [DAS auto zero]                     SL: 0
(0043, 104a) [Starting channel of view]          SS: 1
(0043, 104b) [DAS xm pattern]                    SL: 0
(0043, 104c) [TGGC trigger mode]                 SS: 0
(0043, 104d) [Start scan to X-ray on delay]      FL: 0.0
(0043, 104e) [Duration of X-ray on]              FL: 10.60060977935791
(7fe0, 0010) Pixel Data                          OW: Array of 32768 elements
(fffc, fffc) Data Set Trailing Padding           OB: Array of 126 elements

Plot Images¶

Matplotlib's imshow() function </b> cmap='gray' means gray color mappings for each value

In [7]:
import matplotlib.pyplot as plt
In [8]:
# Draw the image in grayscale
plt.imshow(im_ct, cmap='gray')

# Render the image
plt.show()

Load Volumes¶

ImageIO's volread() function
Load multi-dimensional datasets and create 3D volumes from a folder of images

In [9]:
# Path to the folder containing the DICOM files
folder_path = 'Data/Brain/SE000001/MR000000'

# Load the brain data folder
vol = imageio.volread(folder_path)
Reading DICOM (examining files): 1/27 files (3.7%27/27 files (100.0%)
  Found 1 correct series.
Reading DICOM (loading data): 27/27  (100.0%)

Slice 3D Images¶


Slicing 3D and 4D images into many 2D frames

In [10]:
# Plot the first slice of the volume (2D frame)
plt.imshow(vol[0, :, :], cmap='gray') 

# Render the image
plt.show()

Plot Other Views¶


Slicing along different axes

In [12]:
# Compute aspect ratios
d0, d1, d2 = vol.meta['sampling']
asp1 = d0 / d2
# Plot the other view slice of the volume
plt.imshow(vol[:, 128, :], cmap='gray',aspect=asp1)
# Render the image
plt.show()

Picture2.png

Image Comparison¶

Translation¶


Shifting or moving an image from one place to another within a picture.

In [14]:
import scipy.ndimage as ndi
In [16]:
#formatting method for the plots in this notebook:
def format_and_render_plot():
    '''
    Custom function to simplify common formatting operations for exercises. Operations include: 
    1. Turning off axis grids.
    2. Calling `plt.tight_layout` to improve subplot spacing.
    3. Calling `plt.show()` to render plot.
    '''
    fig = plt.gcf()
    for ax in fig.axes:
        ax.axis('off')    
    plt.tight_layout()
    plt.show()
In [22]:
# save the middle slice as separat image
middle_slice = vol.shape[0] // 2            # // is floor division
im = vol[middle_slice,:,:]

# Translate the brain towards the center
xfm = ndi.shift(im, shift=(10, 15))

# Plot the original and adjusted images
fig, axes = plt.subplots(nrows=1, ncols=2)
axes[0].imshow(im, cmap='gray')
axes[0].set_title('original', fontweight ="bold")
axes[1].imshow(xfm, cmap='gray')
axes[1].set_title('translated', fontweight ="bold")
format_and_render_plot()

Rotation¶


Rotating an image from its center by the specified degrees from the right horizontal axis.

In [18]:
# Rotate the shifted image
xfm = ndi.rotate(xfm, angle=-30, reshape=False)

# Plot the original and transformed images
fig, axes = plt.subplots(nrows=1, ncols=2)
axes[0].imshow(im, cmap='gray')
axes[0].set_title('original', fontweight ="bold")
axes[1].imshow(xfm, cmap='gray')
axes[1].set_title('rotated', fontweight ="bold")
format_and_render_plot()

Affine Transformation¶


Scaling (changing size), rotation, translation (shifting), and shearing (changing angles) using matrix multiplication

affine_transform.png

In [21]:
# Define the affine transform matrix
mat = np.array([[0.8, -0.4, 90], [0.4, 0.8, -6.0], [0, 0, 1]])
print(mat)

# Apply the affine transform matrix to image data
xfm = ndi.affine_transform(im, mat)

# Plot the original and transformed images
fig, axes = plt.subplots(nrows=1, ncols=2)
axes[0].imshow(im, cmap='gray')
axes[0].set_title('original', fontweight ="bold")
axes[1].imshow(xfm, cmap='gray')
axes[1].set_title('transformed', fontweight ="bold")
format_and_render_plot()
[[ 0.8 -0.4 90. ]
 [ 0.4  0.8 -6. ]
 [ 0.   0.   1. ]]

Resampling¶

  • Downsampling: combining pixel data to decrease size
  • Upsampling: distributing pixel data to increase size
In [23]:
# Resample image
im_dn = ndi.zoom(im, zoom=0.25)
im_up = ndi.zoom(im, zoom=4.00)

# Plot the images
fig, axes = plt.subplots(1, 2)
axes[0].imshow(im_dn, cmap='gray')
axes[0].set_title('downsampled', fontweight ="bold")
axes[1].imshow(im_up, cmap='gray')
axes[1].set_title('upsampled', fontweight ="bold")
format_and_render_plot()

Interpolation¶

Estimating new pixel intensities after an image transformation

In [24]:
# Upsample "im" by a factor of 4
up0 = ndi.zoom(im, zoom=512/128, order=0)
up5 = ndi.zoom(im, zoom=512/128, order=5)

# Print original and new shape
print('Original shape:', im.shape)
print('Upsampled shape:', up5.shape)

# Plot close-ups of the new images
fig, axes = plt.subplots(1, 2)
axes[0].imshow(up0[128:256, 128:256], cmap='gray')
axes[0].set_title('no interpolation', fontweight ="bold")
axes[1].imshow(up5[128:256, 128:256], cmap='gray')
axes[1].set_title('with interpolation', fontweight ="bold")
format_and_render_plot()
Original shape: (256, 256)
Upsampled shape: (1024, 1024)

Mean Absolute Error (MAE)¶


Mean intensity differences between two images, with higher values indicating greater divergence

In [25]:
# Load the dcm file (image)
im1 = im

# Apply the affine transform matrix to image data
xfm = ndi.shift(im, shift=(-12, -8))
im2 = xfm

# Calculate image difference
err = im1 - im2

# Plot the difference
fig, axes = plt.subplots(1, 3)
axes[0].imshow(im1, cmap='gray')
axes[0].set_title('original', fontweight ="bold")
axes[1].imshow(im2, cmap='gray')
axes[1].set_title('shifted', fontweight ="bold")
axes[2].imshow(err, cmap='seismic', vmin=-200, vmax=200)
axes[2].set_title('error map', fontweight ="bold")
format_and_render_plot()

# Calculate absolute image difference
abs_err = np.absolute(im1 - im2)

# Plot the difference
fig, axes = plt.subplots(1, 3)
axes[0].imshow(im1, cmap='gray')
axes[0].set_title('original', fontweight ="bold")
axes[1].imshow(im2, cmap='gray')
axes[1].set_title('shifted', fontweight ="bold")
axes[2].imshow(abs_err, cmap='seismic', vmin=-200, vmax=200)
axes[2].set_title('error map (absolute)', fontweight ="bold")
format_and_render_plot()

# Calculate mean absolute error
mean_abs_err = np.mean(np.abs(im1 - im2))
print('MAE:', mean_abs_err)
MAE: 25730.311309814453

Structural Similarity (SSIM)¶

Structural similarity between two images perceived changes in:
structural information, luminance, and contrast that humans typically notice.

In [26]:
from skimage.metrics import structural_similarity as ssim

# Compute SSIM between im and xfm
data_range=xfm.max() - xfm.min()
ssim_index = ssim(im, xfm, data_range=data_range)

print(f'SSIM: {ssim_index}')
SSIM: 0.47951331548119075

Intersection of The Union (IoU)¶

The number of pixels filled in both images (the intersection) out of the number of pixels filled in either image (the union).

In [27]:
def intersection_of_union(im1, im2):
    i = np.logical_and(im1, im2)
    u = np.logical_or(im1, im2)
    return i.sum() / u.sum()

# Try some other paramters by yourself
xfm = ndi.shift(im, shift=(-10, -10))
xfm = ndi.rotate(xfm, angle=-15, reshape=False)
intersection_of_union(xfm, im)
Out[27]:
0.8135092272202998

Combining Two Images¶

The number of pixels filled in both images (the intersection) out of the number of pixels filled in either image (the union).

In [35]:
import SimpleITK as sitk
import Utilities.gui
img1 = sitk.ReadImage('../ImageAnalysis/Data/Brain/training_001_mr_T1.mha') #MRI image
img2_original = sitk.ReadImage('../ImageAnalysis/Data/Brain/training_001_ct.mha')  #CT image
img2 = sitk.Resample(img2_original, img1)

# Obtain foreground masks for the two images using Otsu thresholding, we use these later on.
msk1 = sitk.OtsuThreshold(img1,0,1)
msk2 = sitk.OtsuThreshold(img2,0,1)

Utilities.gui.MultiImageDisplay(image_list = [img1, img2],                   
                      title_list = ['image1', 'image2'],
                      figure_size=(9,3))
Out[35]:
<Utilities.gui.MultiImageDisplay at 0x2a5c2f58280>

Picture3.png

Masks and Filters¶

Load Images¶

In [42]:
#  now we load test data files (CT and MR)
ct_file = get_testdata_file("CT_small.dcm")

# Load the CT file (image)
im = imageio.imread(ct_file)
C:\Users\nasta\AppData\Local\Temp/ipykernel_25900/2938429760.py:5: DeprecationWarning: Starting with ImageIO v3 the behavior of this function will switch to that of iio.v3.imread. To keep the current behavior (and make this warning disappear) use `import imageio.v2 as imageio` or call `imageio.v2.imread` directly.
  im = imageio.imread(ct_file)

Intensity¶

Intensity Normalization¶

In [43]:
# min-max normalisation
im_old = im                                     # save original image

print('Min. value before normalization:', im_old.min())
print('Max value before normalization:', im_old.max())

im = (im - im.min()) / (im.max() - im.min())    # normalise to 0-1 range

print('Min. value after normalization:', im.min())
print('Max value after normalization:', im.max())
Min. value before normalization: -896
Max value before normalization: 1167
Min. value after normalization: 0.0
Max value after normalization: 1.0
In [44]:
def format_and_render_plot(axis=False, legend=False):
    '''
    Custom function to simplify common formatting operations for exercises. Operations include: 
    1. Turning off axis grids and legends, if not explicitly requested.
    2. Calling `plt.tight_layout` to improve subplot spacing.
    3. Calling `plt.show()` to render plot.
    '''
    fig = plt.gcf()
    for ax in fig.axes:
        if not axis:
            ax.axis('off')
        if legend:  
            ax.legend(loc='center right')  
    plt.tight_layout()
    plt.show()
fig, axes = plt.subplots(1, 2, sharex=True)
axes[0].imshow(im_old, cmap='gray')
axes[0].set_title('without normalization', fontweight ="bold")
axes[1].imshow(im, cmap='gray')
axes[1].set_title('with normalization', fontweight ="bold")
format_and_render_plot()

Histogram¶

The distribution of values in your image by binning each element by its intensity then measuring the size of each bin.

CDF (Cumulative distribution function)¶


The frequency with which a given range of pixel intensities occurs.

In [45]:
# Create a 256-bin histogram, binned at each possible value
hist = ndi.histogram(im, min=im.min(), max=im.max(), bins=256)

# Create a cumulative distribution function
cdf = hist.cumsum() / hist.sum()

# Plot the histogram and CDF
fig, axes = plt.subplots(2, 1, sharex=True)         # sharex=True shares the x-axis between the top and bottom subplot
axes[0].plot(hist, label='Histogram')               # label='Histogram' labels this line as "Histogram" for the legend
axes[1].plot(cdf, label='CDF')                      # label='CDF' labels this line as "CDF" for the legend
format_and_render_plot(axis=True, legend=True)      # axis=True turns on axis grids for the plot; legend=True turns on the legend

Create a Mask¶

Binary arrays for removing or selecting specific parts of an image by applying one or more logical operations to an image.

In [46]:
# Create soft tissue (st) and bone masks
# Try different intervals by yourself
mask_st = (im >= 0.2) & (im < 0.52)     # soft tissue mask
mask_bone = im >= 0.53                  # bone mask

# Plot the skin (0) and bone (1) masks
fig, axes = plt.subplots(1,2)
axes[0].imshow(mask_st, cmap='gray')
axes[0].set_title('soft tissue map', fontweight ="bold")
axes[1].imshow(mask_bone, cmap='gray')
axes[1].set_title('bone map', fontweight ="bold")
format_and_render_plot()

Apply a Mask¶

Although masks are binary, they can be applied to images to filter out pixels where the mask is False.

NumPy's where() function is a flexible way of applying masks. It takes three arguments:

np.where(condition, x, y)

condition, x and y can be either arrays or single values. This allows you to pass through original image values while setting masked values to 0.

In [53]:
# Screen out non-bone pixels from "im"
im_bone = np.where(mask_bone, im, 0)  #extract pixels from im which locate in mask_bone!

# Plot masked image and histogram
plt.imshow(im_bone, cmap='gray', label='Bone pixels')
Out[53]:
<matplotlib.image.AxesImage at 0x2a5c6fd0520>

Tune a Mask¶

Imperfect masks can be tuned through the addition and subtraction of pixels. example methods for accomplishing these ends:

  • binary_dilation: Add pixels along edges
  • binary_closing: Dilate then erode, "filling in" holes
In [55]:
# Create and tune bone mask
mask_dilate = ndi.binary_dilation(mask_bone, iterations=5)
mask_closed = ndi.binary_closing(mask_bone, iterations=5)

# Plot masked images
fig, axes = plt.subplots(1,3)
axes[0].imshow(mask_bone)
axes[0].set_title('original', fontweight ="bold")
axes[1].imshow(mask_dilate)
axes[1].set_title('dilated', fontweight ="bold")
axes[2].imshow(mask_closed)
axes[2].set_title('closed', fontweight ="bold")
format_and_render_plot()

Filter Convolutions¶

Transforming images based on intensity values surrounding a pixel

In [56]:
# Set filter weights
weights = [[0.11, 0.11, 0.11],
           [0.11, 0.11, 0.11], 
           [0.11, 0.11, 0.11]]

# Convolve the image with the filter
im_filt = ndi.convolve(im, weights)

# Plot the images
fig, axes = plt.subplots(1,2)
axes[0].imshow(im, cmap='gray')
axes[0].set_title('original', fontweight ="bold")
axes[1].imshow(im_filt, cmap='gray')
axes[1].set_title('with filter', fontweight ="bold")
format_and_render_plot()

Smoothing¶

Improve the signal-to-noise ratio of an image by blurring out small variations in intensity. The Gaussian filter is excellent for this: it is a circular (or spherical) smoothing kernel that weights nearby pixels higher than distant ones.

In [57]:
# Smooth "im" with Gaussian filters
# Try different sigmas  by yourself
im_s1 = ndi.gaussian_filter(im, sigma=1)
im_s3 = ndi.gaussian_filter(im, sigma=3)

# Draw bone masks of each image
fig, axes = plt.subplots(1,3)
axes[0].imshow(im, cmap='gray')
axes[0].set_title('original', fontweight ="bold")
axes[1].imshow(im_s1, cmap='gray')
axes[1].set_title('Gaussian (sigma=1)', fontweight ="bold")
axes[2].imshow(im_s3, cmap='gray')
axes[2].set_title('Gaussian (sigma=3)', fontweight ="bold")
format_and_render_plot()

Detect Edges¶

Applying filters which pattern is a change in intensity along a plane.

In [58]:
# Set weights to detect vertical edges
weights = [[1, 0, -1],
           [1, 0, -1],
           [1, 0, -1]]

# Convolve "im" with filter weights
edges = ndi.convolve(im, weights)

# Draw the image in color
fig, axes = plt.subplots(1,2)
axes[0].imshow(im, cmap='gray')
axes[0].set_title('original', fontweight ="bold")
axes[1].imshow(edges, cmap='seismic', vmin=-0.8, vmax=0.8)
axes[1].set_title('edge filter', fontweight ="bold")
format_and_render_plot()

Picture4.png

Measurement¶

Load Image/Volume¶

In [60]:
# Load the directory (volume)
folder_path = 'Data/Brain/SE000001/MR000000'

# Load the volume
vol = imageio.volread(folder_path)
print(vol.shape)

# save the middle slice as separat image
middle_slice = vol.shape[0] // 2            # // is floor division
im = vol[middle_slice,:,:]
# min-max normalisation
im_old = im                                     # save original image for later
im = (im - im.min()) / (im.max() - im.min())    # normalise the image, range 0 - 1
Reading DICOM (examining files): 1/27 files (3.7%27/27 files (100.0%)
  Found 1 correct series.
Reading DICOM (loading data): 27/27  (100.0%)
(27, 256, 256)

Segmentation¶

Apply a mask and fill any gaps in the mask

In [61]:
# Smooth intensity values
im_filt = ndi.median_filter(im, size=3)         # size = 3 means 3x3x3 neighbourhood

# Select high-intensity pixels
mask_start = np.where(im_filt > 0.3, 1, 0)      # mask_start is a boolean array
mask = ndi.binary_closing(mask_start)           # fill holes

# Label the objects in "mask"
labels, nlabels = ndi.label(mask)               # labels: each object has a unique number, nlabels: number of objects

# Create a `labels` overlay
overlay = np.where(labels > 0, labels, np.nan)

# Use imshow to plot the overlay
fig, axes = plt.subplots(1, 2)
axes[0].imshow(im, cmap='gray')
axes[0].set_title('original', fontweight ="bold")
axes[1].imshow(im, cmap='gray')                             # show image first
axes[1].imshow(overlay, cmap='rainbow', alpha=0.6)          # show overlay second; alpha controls transparency
axes[1].set_title('with segmentation', fontweight ="bold")
format_and_render_plot()

Select Objects¶

Pick up whole sets of pixels at a time

In [62]:
# Label the image "mask"
labels, nlabels = ndi.label(mask)

# Select brain label
brain_val = 2                                           # 2 is the label for the brain (see plot above)
brain_mask = np.where(labels == brain_val, 1, np.nan)   # create brain mask

# Overlay selected label
fig, axes = plt.subplots(1, 2)
axes[0].imshow(im, cmap='gray')
axes[0].set_title('original', fontweight ="bold")
axes[1].imshow(im, cmap='gray')
axes[1].imshow(brain_mask, cmap='rainbow', alpha=0.6)
axes[1].set_title('with segmentation', fontweight ="bold")
format_and_render_plot()

Extract Objects¶

Crop images so that they only include the object of interest.

In [63]:
brain_mask=brain_mask.astype(np.int64)

# Find bounding box
bboxes =ndi.find_objects(brain_mask)                # returns a list of tuples, each tuple has 3 slices
print('Number of objects:', len(bboxes))            # number of objects found
print('Indices for first box:', bboxes[0])          # print indices for first box

# Crop to index 0
im_brain = im[bboxes[0]]                            # crop the original image to the bounding box

# Plot the cropped image
fig, axes = plt.subplots(1, 2)
axes[0].imshow(im, cmap='gray')
axes[0].set_title('original', fontweight ="bold")
axes[1].imshow(im_brain, cmap='gray')
axes[1].set_title('cropped', fontweight ="bold")
format_and_render_plot()
Number of objects: 1
Indices for first box: (slice(50, 234, None), slice(61, 193, None))

Measure Variance¶

Tailor measurements to specific sets of pixels:

In [64]:
# Variance for all pixels
var_all = ndi.variance(vol, labels=None, index=None)            # labels=None means all pixels, index=None means all objects
print('All pixels:', var_all)

# Variance for labeled pixels
var_labels = ndi.variance(vol, labels, index=None)              # labels=labels means only labeled pixels, index=None means all objects
print('Labeled pixels:', var_labels)

# Variance for each object
var_objects = ndi.variance(vol, labels, index=[1,2])            # labels=labels means only labeled pixels, index=[1,2] means only objects 1 and 2
print('Brain matter:', var_objects[0])
print('Other tissue:', var_objects[1])
All pixels: 15281.189220274951
Labeled pixels: 18188.560688187546
Brain matter: 13466.402775068886
Other tissue: 16319.82340391041

Separate Histograms¶

Wide distribution of intensity values and more variance of multiple tissue types

In [65]:
# Create histograms for selected pixels
hist1 = ndi.histogram(vol, min=0, max=255, bins=256)
hist2 = ndi.histogram(vol, 0, 255, 256, labels=labels)
hist3 = ndi.histogram(vol, 0, 255, 256, labels=labels, index=1)


# Plot the histogram and CDF
fig, axes = plt.subplots(3, 1, sharex=True)                 # sharex=True shares the x-axis between the top and bottom subplot
axes[0].plot(hist1 / hist1.sum(), label='All pixels')
axes[1].plot(hist2 / hist2.sum(), label='All labeled pixels')
axes[2].plot(hist3 / hist3.sum(), label='Brain matter')
format_and_render_plot(axis=True, legend=True)              # axis=True turns on axis grids for the plot; legend=True turns on the legend

Calculate Distance¶

the distance from each pixel to a given point

In [66]:
# Calculate distances
brain = np.where(labels == 2, 1, 0)                                             # create brain mask
dists = ndi.distance_transform_edt(brain, sampling=vol.meta['sampling'][1:3])   # calculate distances

# Report on distances
print('Max distance (mm):', ndi.maximum(dists))
print('Max location:', ndi.maximum_position(dists))

# Plot overlay of distances
overlay = np.where(dists > 0, dists, np.nan)                                    # create overlay

fig, axes = plt.subplots(1, 2)
axes[0].imshow(im, cmap='gray')
axes[0].set_title('original', fontweight ="bold")
axes[1].imshow(im, cmap='gray')
axes[1].imshow(overlay, cmap='hot', alpha=0.75)
axes[1].set_title('with distances', fontweight ="bold")
format_and_render_plot()
Max distance (mm): 12.1875
Max location: (133, 175)

Center-Of-Mass (COM)¶

The coordinates for the point of an object with higher values pulling the center closer to it.

In [67]:
# Extract centers of mass for objects 1 and 2
coms = ndi.center_of_mass(vol, labels, index=[1,2])
print('Label 1 center:', coms[0])
print('Label 2 center:', coms[1])

# Add marks to plot
for c0, c1, c2 in coms:
    plt.scatter(c2, c1, s=100, marker='o')
plt.show()
Label 1 center: (8.856682792208547, 49.44438084201855, 125.03554483621987)
Label 2 center: (11.713445688803832, 119.3292456288401, 129.36906622106693)
In [ ]: